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Abstract. This paper describes a numerical method for finding good packings in Grassmannian 
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manifolds equipped with various metrics. This investigation also encompasses packing in projective 
spaces. In each case, producing a good packing is equivalent to constructing a matrix that has 
certain structural and spectral properties. By alternately enforcing the structural condition and 
Oh. then the spectral condition, it is often possible to reach a matrix that satisfies both. One may then 

' extract a packing from this matrix. 

This approach is both powerful and versatile. In cases where experiments have been performed, 
the alternating projection method yields packings that compete with the best packings recorded. 
It also extends to problems that have not been studied numerically. For example, it can be used to 
produce packings of subspaces in real and complex Grassmannian spaces equipped with the Fubini- 
Study distance; these packings are valuable in wireless communications. One can prove that some 
of the novel configurations constructed by the algorithm have packing diameters that are nearly 
optimal. 
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1. Introduction 



Let us begin with the standard facetious example. Imagine that several mutually inimical nations 
build their capital cities on the surface of a featureless globe. Being concerned about missile strikes, 
I they wish to locate the closest pair of cities as far apart as possible. In other words, what is the 

■ best way to pack points on the surface of a two-dimensional sphere? 

^ , This question, first discussed by the Dutch biologist Tammes |Tam30] , is the prototypical example 

On I of packing in a compact metric space. It has been studied in detail for the last 75 years. More 

■ recently, researchers have started to ask about packings in other compact spaces. In particular, 
I several communities have investigated how to arrange subspaces in a Euclidean space so that they 



are as distinct as possible. An equivalent formulation is to find the best packings of points in 
a Grassmannian manifold. This problem has applications in quantum computing and wireless 
^ I communications. There has been theoretical interest in subspace packing since the 1960s |T6t65j . 

■ but the first detailed numerical study appears in a 1996 paper of Conway, Hardin, and Sloane 

|CHS96j . 

The aim of this paper is to describe a flexible numerical method that can be used to construct 
packings in Grassmannian manifolds equipped with several different metrics. The rest of this 
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introduction provides a formal statement of abstract packing problems, and it offers an overview 
of our approach to solving them. 

1.1. Abstract Packing Problems. Although we will be working with Grassmannian manifolds, 
it is more instructive to introduce packing problems in an abstract setting. Let M be a compact 
metric space endowed with the distance function distM- The packing diameter of a finite subset St^ 
is the minimum distance between some pair of distinct points drawn from ^ . That is, 

packM(^) = min distM (a^mjXn). 

In other words, the packing diameter of a set is the diameter of the largest open ball that can be 
centered at each point of the set without encompassing any other point. (It is also common to 
study the packing radius, which is half the diameter of this ball.) An optimal packing of N points 
is an ensemble that solves the mathematical program 

max packivif^^) 

\X\=N 

where |-| returns the cardinality of a finite set. The optimal packing problem is guaranteed to have 
a solution because the metric space is compact and the objective is a continuous function of the 
ensemble ^ . 

This article focuses on a feasibility problem closely connected with optimal packing. Given a 
number p, the goal is to produce a set of N points for which 

packM(.^) > p. (1.1) 

This problem is notoriously difficult to solve because it is highly nonconvex, and it is even more 
difficult to determine the maximum value of p for which the feasibility problem is soluble. This 
maximum value of p corresponds with the diameter of an optimal packing. 

1.2. Alternating Projection. We will attempt to solve the feasibility problem (jl.ip in Grass- 
mannian manifolds equipped with a number of different metrics, but the same basic algorithm 
applies in each case. Here is a high-level description of our approach. 

First, we show that each configuration of subspaces is associated with a block Gram matrix whose 
blocks control the distances between pairs of subspaces. Then we prove that a configuration solves 
the feasibility problem (jl.ip if and only if its Gram matrix possesses both a structural property 
and a spectral property. The overall algorithm consists of the following steps. 

(1) Choose an initial configuration and construct its matrix. 

(2) Alternately enforce the structural condition and the spectral condition in hope of reaching 
a matrix that satisfies both. 

(3) Extract a configuration of subspaces from the output matrix. 

In our work, we choose the initial configuration randomly and then remove similar subspaces 
from it with a simple algorithm. One can imagine more sophisticated approaches to constructing 
the initial configuration. 

Flexibility and ease of implementation are the major advantages of alternating projection. This 
article demonstrates that appropriate modifications of this basic technique allow us to construct 
solutions to the feasibility problem in Grassmannian manifolds equipped with various metrics. 
Some of these problems have never been studied numerically, and the experiments point toward 
intriguing phenomena that deserve theoretical attention. Moreover, we believe that the possibilities 
of this method have not been exhausted and that it will see other applications in the future. 

Alternating projection does have several drawbacks. It may converge very slowly, and it does 
not always yield a high level of numerical precision. In addition, it may not deliver good packings 
when the ambient dimension or the number of subspaces in the configuration is large. 
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1.3. Motivation and Related Work. This work was motivated by applications in electrical 
engineering. In particular, subspace packings solve certain extremal problems that arise in multiple- 
antenna communication systems |ZT02l iHMR+OOl ILJSH04j . This appHcation requires complex 



Grassmannian packings that consist of a small number of subspaces in an ambient space of low 
dimension. Our algorithm is quite effective in this parameter regime. The resulting packings fill 
a significant gap in the literature, since existing tables consider only the real case |Slo04aj . See 
Section 16.11 for additional discussion of the wireless application. 

The approach to packing via alternating projection was discussed in a previous publication 
[TDJS05] . but the experiments were limited to a single case. We are aware of several other numer- 
ical methods that can be used to construct packings in Grassmannian manifolds |CHS96l ITroOlj 
lARUOlj . These techniques rely on ideas from nonlinear programming. 

1.4. Historical Interlude. The problem of constructing optimal packings in various metric spaces 
has a long and lovely history. The most famous example may be Kepler's Conjecture that an 
optimal packing of spheres in three-dimensional Euclidean spac^ locates them at the points of a 
face-centered cubic lattice. For millennia, greengrocers have applied this theorem when stacking 
oranges, but it has only been established rigorously within the last few years |Hal04| . Packing 
problems play a major role in modern communications because error-correcting codes may be 
interpreted as packings in the Hamming space of binary strings |CT91j . The standard reference on 
packing is the magnum opus of Conway and Sloane |CS98j . Classical monographs on the subject 
were written by L. Fejes Toth |T6t64] and C. A. Rogers |Rog64| . 

The idea of applying alternating projection to feasibility problems first appeared in the work of 
von Neumann |vN50] . He proved that an alternating projection between two closed subspaces of a 
Hilbert space converges to the orthogonal projection of the initial iterate onto the intersection of the 
two subspaces. Cheney and Goldstein subsequently showed that an alternating projection between 
two closed, convex subsets of a Hilbert space always converges to a point in their intersection 
(provided that the intersection is nonempty) |CG59j . This result does not apply in our setting 
because one of the constraint sets we define is not convex. 

1.5. Outline of Article. Here is a brief overview of this article. In Section [2l we develop a 
basic description of Grassmannian manifolds and present some natural metrics. Section [3] explains 
why alternating projection is a natural algorithm for producing Grassmannian packings, and it 
outlines how to apply this algorithm for one specific metric. Section [3] gives some theoretical upper 
bounds on the optimal diameter of packings in Grassmannian manifolds. Section [5] describes the 
outcomes of an extensive set of numerical experiments and explains how to apply the algorithm 
to other metrics. Section [6] offers some discussion and conclusions. Appendix [Al explores how our 
methodology applies to Tammes' Problem of packing on the surface of a sphere. Finally, Appendix 
IB] contains tables and figures that detail the experimental results. 

2. Packing in Grassmannian Manifolds 

This section introduces our notation and a simple description of the Grassmannian manifold. It 
presents several natural metrics on the manifold, and it shows how to represent a configuration of 
subspaces in matrix form. 

2.1. Preliminaries. We work in the vector space C^. The symbol * denotes the complex-conjugate 
transpose of a vector (or matrix). We equip the vector space with its usual inner product (a;, y) = 
y*x. This inner product generates the £2 norm via the formula ||a;||2 = {x, x). 

The d-dimensional identity matrix is 1^; we sometimes omit the subscript if it is unnecessary. A 
square matrix is positive semidefinite when its eigenvalues are all nonnegative. We write X ;^ to 
indicate that X is positive semidefinite. 



-'^The infinite extent of a Euclidean space necessitates a more subtle definition of an optimal packing. 
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A square, complex matrix U is unitary if it satisfies U*U = I. If in addition the entries of U 
are real, the matrix is orthogonal. The unitary group U(d) can be presented as the collection of all 
dx d unitary matrices with ordinary matrix multiplication. The real orthogonal group 0{d) can be 
presented as the collection of all d x d real orthogonal matrices with the usual matrix multiplication. 

Suppose that X is a general matrix. The Frobenius norm is calculated as ||-X^||f = trace X* X , 
where the trace operator sums the diagonal entries of the matrix. The spectral norm is denoted by 
II II2 2; it returns the largest singular value of X. Both these norms are unitarily invariant, which 
means that ||C/XV*|| = ||^|| whenever U and V are unitary. 

2.2. Grassmannian Manifolds. The (complex) Grassmannian manifold G{K, C^) is the collec- 
tion of all A'-dimensional subspaces of C^. This space is isomorphic to a quotient of unitary groups: 

U(d) 



G{K, 



U{K) X U((i - K) ' 



To understand the equivalence, note that each orthonormal basis from can be split into K 
vectors, which span a i^T-dimensional subspace, and d — K vectors, which span the orthogonal 
complement of that subspace. To obtain a unique representation for the subspace, it is necessary 
to divide by isometrics that fix the subspace and by isometrics that fix its complement. It is evident 
that G{K,C^) is always isomorphic to G{d — K,C^). 

Similarly, the real Grassmannian manifold G{K, W^) is the collection of all i^-dimensional sub- 
spaces of M*^. This space is isomorphic to a quotient of orthogonal groups: 



0{K) X 0{d - K) ' 

If we need to refer to the real and complex Grassmannians simultaneously, we write G{K,¥'^). 

In the theoretical development, we concentrate on complex Grassmannians since the development 
for the real case is identical, except that all the matrices are real-valued instead of complex- valued. 
A second reason for focusing on the complex case is that complex packings arise naturally in wireless 
communications |LJS03j . 

When each subspace has dimension = 1, the Grassmannian manifold reduces to a simpler 
object called a projective space. The elements of a projective space can be viewed as lines through 
the origin of a Euclidean space. The standard notation is P'^~'^(F) G(1,F'^). We will spend a 
significant amount of attention on packings of this manifold. 

2.3. Principal Angles. Suppose that S and T are two subspaces in G{K,G'^). These subspaces 
are inclined against each other by K different principal angles. The smallest principal angle 9i is 
the minimum angle formed by a pair of unit vectors (si, ti) drawn from S x T. That is, 

9i = min arccos(si, ti) subject to ||si|L = 1 and ||ti|L = 1. 
(si,ti)e5xr 

The second principal angle 92 is defined as the smallest angle attained by a pair of unit vectors 
(^2,^2) that is orthogonal to the first pair, i.e., 

92= min arccos(s2; ^2) subject to IIS2II2 = 1 P2II2 = 

{s2,t2)eSxT 

{si, 82)= and (ii, i2) = 0. 

The remaining principal angles are defined analogously. The sequence of principal angles is nonde- 
creasing, and it is contained in the range [0,7r/2]. We only consider metrics that are functions of 
the principal angles between two subspaces. 

Let us present a more computational definition of the principal angles |BG73j . Suppose that the 
columns of S and T form orthonormal bases for the subspaces S and T. More rigorously, is a 
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dx K matrix that satisfies S*S = Ik and range S = S. The matrix T has an analogous definition. 
Next we compute a singular value decomposition of the product S*T: 

S*T = UCV* 

where U and V are K x K unitary matrices and C is a nonnegative, diagonal matrix with non- 
increasing entries. The matrix C of singular values is uniquely determined, and its entries are the 
cosines of the principal angles between S and T: 

Ckk = cos6k k = l,2,...,K. 

This definition of the principal angles is most convenient numerically because singular value de- 
compositions can be computed efficiently with standard software. We also note that this definition 
of the principal angles does not depend on the choice of matrices S and T that represent the two 
subspaces. 

2.4. Metrics on Grassmannian Manifolds. Grassmannian manifolds admit many interesting 
metrics, which lead to different packing problems. This section describes some of these metrics. 

(1) The chordal distance between two i^T-dimensional subspaces S and T is given by 



distchord('5,T) = V sin^ 6*1 H h sin^ Ok 



K -\\S*T\\l 



1/2 

(2.1) 



The values of this metric range between zero and ^/K. The chordal distance is the easiest 
to work with, and it also yields the most symmetric packings |CHS96| . 
(2) The spectral distance is 



distspec(5, T) = miufc sin 0k 



\S*T\\l, 



1/2 

(2.2) 



The values of this metric range between zero and one. As we will see, this metric promotes 
a special type of packing called an equi-isoclinic configuration of subspaces. 
(3) The Fubini-Study distance is 



distFs('5, T) = arccos (^JJ^ cos ( 



= arccos|detS'*T| . (2.3) 

This metric takes values between zero and 7r/2. It plays an important role in wireless 
communications |LHJ051 ILJOSj . 

(4) The geodesic distance is 



distgeo(5,r)'=Ve^ + --- + 0|. 



This metric takes values between zero and 7ry/K/2. From the point of view of differen- 
tial geometry, the geodesic distance is very natural, but it does not seem to lead to very 
interesting packings |CHS96j . so we will not discuss it any further. 

Grassmannian manifolds support several other interesting metrics, some of which are listed in 
[BN02j . In case we are working in a projective space, i.e., K = 1, all of these metrics reduce to 
the acute angle between two lines or the sine thereof. Therefore, the metrics are equivalent up to 
a monotonically increasing transformation, and they promote identical packings. 
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2.5. Representing Configurations of Subspaces. Suppose that ^ = {Si, . . . ,Sn} is a col- 
lection of N subspaces in G(i^, C^). Let us develop a method for representing this configuration 
numerically. To each subspace we associate a (nonunique) d x K matrix Xn whose columns 
form an orthonormal basis for that subspace, i.e., X*Xn = Ix and range X„ = Now collate 
these matrices into a d x KN configuration matrix 

X^=*[Xi X.2 ... Xn]. 

In the sequel, we do not distinguish between the configuration and the matrix X. 

The Gram matrix of X is defined as the KN x KN matrix G = X*X. By construction, the 
Gram matrix is positive semidefinite, and its rank does not exceed d. It is best to regard the Gram 
matrix as an x block matrix comprised oi K x K blocks, and we index it as such. Observe 
that each block satisfies 

In particular, each diagonal block Gnn is an identity matrix. Meanwhile, the singular values of 
the off-diagonal block Gmn equal the cosines of the principal angles between the two subspaces 
range and range 

Conversely, let G be an x block matrix with each block of size K x K. Suppose that the 
matrix is positive semidefinite, that its rank does not exceed d, and that its diagonal blocks are 
identity matrices. Then we can factor G = X*X where X is a dx KN configuration matrix. That 
is, the columns of X form orthogonal bases for A^ different X-dimensional subspaces of C"'. 

As we will see, each metric on the Grassmannian manifold leads to a measure of "magnitude" for 
the off-diagonal blocks on the Gram matrix G. A configuration solves the feasibility problem (II. ip 
if and only if each off-diagonal block of its Gram matrix has sufficiently small magnitude. So solving 
the feasibility problem is equivalent to producing a Gram matrix with appropriate properties. 

3. Alternating Projection for Chordal Distance 

In this section, we elaborate on the idea that solving the feasibility problem is equivalent with 
constructing a Gram matrix that meets certain conditions. These conditions fall into two different 
categories: structural properties and spectral properties. This observation leads naturally to an 
alternating projection algorithm for solving the feasibility problem. The algorithm alternately 
enforces the structural properties and then the spectral properties in hope of producing a Gram 
matrix that satisfies them all. This section illustrates how this approach unfolds when distances 
are measured with respect to the chordal metric. In Section [SJ we describe adaptations for other 
metrics. 

3.1. Packings with Chordal Distance. Suppose that we seek a packing of A^ subspaces in 
G{K, C^) equipped with the chordal distance. If X is a configuration of A^ subspaces, its packing 
diameter is 



packchoj.d(^) = min distchord (^m,^r, 



mm 



^ ~ ll^m^nllF 



1/2 



Given a parameter p, the feasibility problem elicits a configuration X that satisfies 



mm 



^ ~ ll^m^nllp 



1/2 

>P- 



We may rearrange this inequality to obtain a simpler condition: 



max ||X*„X„||p < /i (3.1) 
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where 

H = - p2. (3.2) 

In fact, we may formulate the feasibihty problem pm'ely in terms of the Gram matrix. Suppose 
that the configuration X satisfies (j3.ip with parameter Then its Gram matrix G must have the 
following six properties: 

(1) G is Hermitian. 

(2) Each diagonal block of G is an identity matrix. 

(3) ||Gmn.|lF — M each m ^ n. 

(4) G is positive semidefinite. 

(5) G has rank d or less. 

(6) G has trace KN. 

Some of these properties are redundant, but we have listed them separately for reasons soon to 
become apparent. Conversely, suppose that a matrix G satisfies Properties 1-6. Then it is always 
possible to factor it to extract a configuration of N subspaces that solves (13. ip . The factorization 
of G = X*X can be obtained most easily from an eigenvalue decomposition of G. 

3.2. The Algorithm. Observe that Properties 1-3 are structural properties. By this, we mean 
that they constrain the entries of the Gram matrix directly. Properties 4-6, on the other hand, are 
spectral properties. That is, they control the eigenvalues of the matrix. It is not easy to enforce 
structural and spectral properties simultaneously, so we must resort to half measures. Starting 
from an initial matrix, our algorithm will alternately enforce Properties 1-3 and then Properties 
4-6 in hope of reaching a matrix that satisfies all six properties at once. 
To be more rigorous, let us define the structural constraint set 

J^ifi) = {H e C^^x^^ : H = H*,Hnn = lK iov n = l,2,...,N, 

and llH'mnllp < fi for all m ^ n}. (3.3) 

Although the structural constraint set evidently depends on the parameter fi, we will usually 
eliminate fi from the notation for simplicity. We also define the spectral constraint set 

^ = {G G C-^^''^^ : G ^ 0,rankG < d, and trace G = KN} . (3.4) 

Both constraint sets are closed and bounded, hence compact. The structural constraint set is 
convex, but the spectral constraint set is not. 

To solve the feasibility problem (13. ip , we must find a matrix that lies in the intersection of ^ and 
J^. This section states the algorithm, and the succeeding two sections provide some implementation 
details. 

Algorithm 1 (Alternating Projection). 
Input: 

• A KN X KN Hermitian matrix G^^^ 

• The maximum number of iterations T 

Output: 

• A KN X KN matrix Gout that belongs to and whose diagonal blocks are identity matrices 
Procedure: 

(1) Initialize t <— 0. 

(2) Determine a matrix H^^^ that solves 

min ||i3"-GW|L. 
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(3) Determine a matrix G^*"*"^) that solves 

mm||G-i3"W|L. 

(4) Increment t. 

(5) Ift<T, return to Step 2. 

(6) Define the block-diagonal matrix D = diagG^^^. 

(7) Return the matrix 

Gout = £>-l/2g,(T)^-l/2_ 

The iterates generated by this algorithm are not guaranteed to converge in norm. Therefore, we 
have chosen to halt the algorithm after a fixed number of steps instead of checking the behavior of 
the sequence of iterates. We discuss the convergence properties of the algorithm in the sequel. 

The scaling in the last step normalizes the diagonal blocks of the matrix but preserves its inertia 
(i.e., numbers of negative, zero, and positive eigenvalues). Since is a positive-semidefinite 

matrix with rank d or less, the output matrix Gout shares these traits. It follows that the output 
matrix always admits a factorization Gout = X*X where X is a d x KN configuration matrix. 
Property 3 is the only one of the six properties that may be violated. 

3.3. The Matrix Nearness Problems. To implement Algorithm [H we must solve the matrix 
nearness problems in Steps 2 and 3. The first one is straightforward. 

Proposition 2. Let G he an Hermitian matrix. With respect to the Frohenius norm, the unique 
matrix in ^(/x) nearest to G has diagonal blocks equal to the identity and off-diagonal blocks that 
satisfy 

H ={ ^""^ if\\Gmn\\Y<l^'0.'^d 

\ ^lGmn/ \\Gmn\\Y Otherwise. 

It is rather more difficult to find a nearest matrix in the spectral constraint set. To state the 
result, we define the plus operator by the rule = max{0,a;}. 

Proposition 3. Let H be an Hermitian matrix whose eigenvalue decomposition is Ylf^ '^j'^j'^j 
with the eigenvalues arranged in nonincreasing order: Ai > A2 > • • • > Aa'a^- With respect to the 
Frobenius norm, a matrix in ^ closest to H is given by 

where the scalar 7 is chosen so that 

This best approximation is unique provided that Xa > A^+i . 

The nearest matrix described by this theorem can be computed efficiently from an eigenvalue 
decomposition of H. (See |GVL96j for computational details.) The value of 7 is uniquely deter- 
mined, but one must solve a small rootfinding problem to solve it. The bisection method is an 
appropriate technique since the plus operator is nondifferentiable. We omit the details, which are 
routine. 

Proof. Given an Hermitian matrix A, denote by A(A) the vector of eigenvalues arranged in nonin- 
creasing order. Then we may decompose A = U {diag X{A)}U* for some unitary matrix U. 
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Finding the matrix in closest to H. is equivalent to solving the optimization problem 

min ||G — ||p subject to ^j{G) > for j = 1, . . . ,d, 
G 

Xj{G) = for i = d + 1, . . . , KN, and 



First, we fix the eigenvalues of G and minimize with respect to the unitary part of its eigenvalue 
decomposition. In consequence of the Hoffman- Wielandt Theorem |HJ85| . the objective function 
is bounded below: 



Equality holds if and only if G and H are simultaneously diagonalizable by a unitary matrix. 
Therefore, if we decompose H = U {diag X{H)}U* , the objective function attains its minimal 
value whenever G = U {diag X{G)}U* . Note that the matrix U may not be uniquely determined. 

We find the optimal vector of eigenvalues ^ for the matrix G by solving the (strictly) convex 
program 



This minimization is accomplished by an application of Karush-Kuhn-Tucker theory [Roc70j . In 
short, the top d eigenvalues of H are translated an equal amount, and those that become negative 
are set to zero. The size of the translation is chosen to fulfill the third condition (which controls 
the trace of G). The entries of the optimal ^ are nonincreasing on account of the ordering of X{H). 

Finally, the uniqueness claim follows from the fact that the eigenspace associated with the top d 
eigenvectors of H is uniquely determined if and only if \d[H) > X(i+i{H). □ 

3.4. Choosing an Initial Configuration. The success of the algorithm depends on adequate 
selection of the input matrix G^'^\ We have found that the following strategy is reasonably effective. 
It chooses random subspaces and adds them to the initial configuration only if they are sufficiently 
distant from the subspaces that have already been chosen. 

Algorithm 4 (Initial Configuration). 
Input: 

• The ambient dimension d, the subspace dimension K , and the number N of subspaces 

• An upper bound r on the similarity between subspaces 

• The maximum number T of random selections 

Output: 

• A KN X KN matrix G from 'S whose off-diagonal blocks also satisfy \\Gmn\\F ^ 7" 
Procedure: 

(1) Initialize t <— and n ^ 1. 

(2) Increment t. If t > T, print a failure notice and stop. 

(3) Pick a d X K matrix whose range is a uniformly random subspace in G{K,C^). 

(4) If \\X^Xn\\p < T for each m = 1, . . . ,n — 1, then increment n. 

(5) If n < N , return to Step 2. 

(6) Form the matrix X = \_Xi X2 ■ ■ ■ X]\f~^ . 

(7) Return the Gram matrix G = X* X . 



G-H\\l>\\X(G)-\{H)\\l. 



nun 11^ — A(i3")||2 subject to > for j = 1, . . . , d, 

= for j = d + 1, . . . , KN, and 
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To implement Step 3, we use the method developed in |Ste80j . Draw a d x K matrix whose 
entries are iid complex, standard normal random variables, and perform a QR decomposition. The 
first K columns of the unitary part of the QR decomposition form an orthonormal basis for a 
random X-dimensional subspace. 

The purpose of the parameter r is to prevent the starting configuration X from containing blocks 
that are nearly identical. The extreme case r = ^/K places no restriction on the similarity between 
blocks. If r is chosen too small (or if we are unlucky in our random choices), then this selection 
procedure may fail. For this reason, we add an iteration counter to prevent the algorithm from 
entering an infinite loop. We typically choose values of r very close to the maximum value. 

3.5. Theoretical Behavior of Algorithm. It is important to be aware that packing problems 
are typically difficult to solve. Therefore, we cannot expect that our algorithm will necessarily 
produce a point in the intersection of the constraint sets. One may ask whether we can make 
any guarantees about the behavior of Algorithm [TJ This turns out to be difficult. Indeed, there 
is potential that an alternating projection algorithm will fail to generate a convergent sequence 
of iterates |Mey76| . Nevertheless, it can be shown that the sequence of iterates has accumulation 
points and that these accumulation points satisfy a weak structural property. 

In practice, the alternating projection algorithm seems to converge, but a theoretical justification 
for this observation is lacking. A more serious problem is that the algorithm frequently requires as 
many as 5000 iterations before the iterates settle down. This is one of the major weaknesses of our 
approach. 

For reference, we offer the best theoretical convergence result that we know. The distance 
between a matrix and a compact collection of matrices is defined as 

dist(M,'^) = min \\M - C\\-p . 
It can be shown that the distance function is Lipschitz, hence continuous. 

Theorem 5 (Global Convergence). Suppose that Algorithm d generates an infinite sequence of 
iterates {(GW,//^)}. This sequence has at least one accumulation point. 

• Every accumulation point lies in^ x M' . 

• Every accumulation point (G, H) satisfies 

\\G-H\L= lim IIgW -/i"(*)|L. 

• Every accumulation point (G, H) satisfies 

\\G-H\\^ = dist(G,^) = dist(S^,^). 

Proof sketch. The existence of an accumulation point follows from the compactness of the constraint 
sets. The algorithm does not increase the distance between successive iterates, which is bounded 
below by zero. Therefore, this distance must converge. The distance functions are continuous, so 
we can take limits to obtain the remaining assertions. □ 

A more detailed treatment requires the machinery of point-to-set maps, and it would not enhance 
our main discussion. Please see the appendices of |TDJS05"] for additional information. 

4. Bounds on the Packing diameter 

To assay the quality of the packings that we produce, it helps to have some upper bounds on 
the packing diameter. If a configuration of subspaces has a packing diameter close to the upper 
bound, that configuration must be a nearly optimal packing. This approach allows us to establish 
that many of the packings we construct numerically have packing diameters that are essentially 
optimal. 
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Theorem 6 (Conway-Hardin-Sloane |CHS96j ). The packing diameter of N subspaces in the Grass- 
mannian manifold G(Er, F'^) equipped with chordal distance is bounded above as 

,^,0 K(d-K) N 
pacWd(=^)' < -J^- (4.1) 

If the bound is met, all pairs of subspaces are equidistant. When F = M, the bound is attainable 
only if N < + 1). When F = C, the bound is attainable only if N < d"^. 

The complex case is not stated in |CHS96j . but it fohows from an identical argument. We refer to 
(j4.ip as the Rankin bound for subspace packings with respect to the chordal distance. The reason for 
the nomenclature is that the result is established by embedding the chordal Grassmannian manifold 
into a Euclidean sphere and applying the classical Rankin bound for sphere packing |Ran47] . 

It is also possible to draw a corollary on packing with respect to the spectral distance; this result 
is novel. A subspace packing is said to be equi-isoclinic if all the principal angles between all pairs 
of subspaces are identical |LS73] . 

Corollary 7. We have the following bound on the packing diameter of N subspaces in the Grass- 
mannian manifold G{K, W^) equipped with the spectral distance. 



pack3p,,(jr)2 < 'L^JL_, (4.2) 



If the bound is met, the packing is equi-isoclinic. 

We refer to (j4.2p as the Rankin bound for subspace packings with respect to spectral distance. 

Proof. The power mean inequality (equivalently, Holder's inequality) yields 



miufc sin 0^ < 



1 1/2 

K-iy" sin^^fc , 



For angles between zero and vr/2, equality holds if and only \i 9i = ■ ■ ■ = Ok- It follows that 

pack (jr)2 < Jf-^pack^hordC-^)^ ^d K N 



d N -I 



If the second inequality is met, then all pairs of subspaces are equidistant with respect to the 
chordal metric. Moreover, if the first inequality is met, then the principal angles between each 
pair of subspaces are constant. Together, these two conditions imply that the packing is equi- 
isoclinic. □ 

An upper bound on the maximum number of equi-isoclinic subspaces is available. Its authors do 
not believe that it is sharp. 

Theorem 8 (Lemmens-Seidel |LS73] ). The maximum number of equi-isoclinic K-dimensional 
subspaces ofW^ is no greater than 

\d{d+l)-\K{K + l) + l. 
Similarly, the maximum number of equi-isoclinic K-dimensional subspaces of C"' does not exceed 

d^ -K"^ + 1. 



CONSTRUCTING GRASSMANNIAN PACKINGS 



12 



5. Experiments 

Our approach to packing is experimental rather than theoretical, so the real question is how 
Algorithm [1] performs in practice. In principle, this question is difficult to resolve because the 
optimal packing diameter is unknown for almost all combinations of d and N. Whenever possible, 
we compared our results with the Rankin bound and with the "world record" packings tabulated 
by N. J. A. Sloane and his colleagues |Slo04aj . In many cases, the algorithm was able to identify a 
nearly optimal packing. Moreover, it yields interesting results for packing problems that have not 
received numerical attention. 

In the next subsection, we describe detailed experiments on packing in real and complex pro- 
jective spaces. Then, we move on to packings of subspaces with respect to the chordal distance. 
Afterward, we study the spectral distance and the Fubini-Study distance. 

5.1. Projective Packings. Line packings are the simplest type of Grassmannian packing, so they 
offer a natural starting point. Our goal is to produce the best packing of lines in P'^~^(F). In 
the real case, Sloane's tables allow us to determine how much our packings fall short of the world 
record. In the complex setting, there is no comparable resource, so we must rely on the Rankin 
bound to gauge how well the algorithm performs. 

Let us begin with packing in real projective spaces. We attempted to construct configurations of 
real lines whose maximum absolute inner product /x fell within 10~^ of the best value tabulated in 
|Slo04a] . For pairs {d, N) with (i = 3, 4, 5 and A" = 4, 5, . . . , 25, we computed the putatively optimal 
value of the feasibility parameter fi from Sloane's data and equation (|3.2p . In each of 10 trials, 
we constructed a starting matrix using Algorithm d] with parameters r = 0.9 and T = 10, 000. 
(Recall that the value of T determines the maximum number of random subspaces that are drawn 
when trying to construct the initial configuration.) We applied alternating projection. Algorithm 
[H with the computed value of n and the maximum number of iterations T = 5000. (Our numerical 
experience indicates that increasing the maximum number of iterations beyond 5000 does not confer 
a significant benefit.) We halted the iteration in Step 4 if the iterate G^*^ exhibited no off-diagonal 
entry with absolute value greater than + 10~^. After 10 trials, we recorded the largest packing 
diameter attained, as well as the average value of the packing diameter. We also recorded the 
average number of iterations the alternating projection required per trial. 

Table [T] delivers the results of this experiment. Following Sloane, we have reported the degrees 
of arc subtended by the closest pair of lines. We believe that it is easiest to interpret the results 
geometrically when they are stated in this fashion. All the tables and figures related to packing are 
collated at the back of this paper for easy comparison. 

According to the table, the best configurations produced by alternating projection consistently 
attain packing diameters tenths or hundredths of a degree away from the best configurations known. 
The average configurations returned by alternating projection are slightly worse, but they usually 
fall within a degree of the putative optimal. Moreover, the algorithm finds certain configurations 
with ease. For the pair (5,16), fewer than 1000 iterations are required on average to achieve a 
packing within 0.001 degrees of optimal. 

A second observation is that the alternating projection algorithm typically performs better when 
the number N of points is small. The largest errors are all clustered at larger values of N. A 
corollary observation is that the average number of iterations per trial tends to increase with the 
number of points. 

There are several anomalies that we would like to point out. The most interesting pathology 
occurs at the pair {d, N) = (5, 19). The best packing diameter calculated by alternating projection 
is about 1.76° worse than the optimal configuration, and it is also 1.76° worse than the best packing 
diameter computed for the pair (5, 20). From Sloane's tables, we can see that the (putative) optimal 
packing of 19 lines in P^(]R) is actually a subset of the best packing of 20 lines. Perhaps the fact that 
this packing is degenerate makes it difficult to construct. A similar event occurs (less dramatically) 
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at the pair (5,13). The table also shows that the algorithm performs less effectively when the 
number of lines exceeds 20. 

In complex projective spaces, this methodology does not apply because there are no tables 
available. In fact, we only know of one paper that contains numerical work on packing in complex 
projective spaces [ARUOlj . but it gives very few examples of good packings. The only method we 
know for gauging the quality of a complex line packing is to compare it against an upper bound. The 
Rankin bound for projective packings, which is derived in Section [H states that every configuration 
^ of hues in either P'^-i(]R) or P'^-i(C) satisfies the inequality 

oack i3t:? < (^-^^^ 
packp(^) <^(^_-^)- 

This bound is attainable only for rare combinations of d and N. In particular, the bound can be 
met in P'^~^(M) only if TV < i d (d + 1). In the space P'='-i(C), attainment requires that N < (P. 
Any arrangement of lines that meets the Rankin bound must be equiangular. These optimal 
configurations are called equiangular tight frames. See |SJ031 IIIP041 ITDJSOSl ISTDJ07] for more 
details. 

We performed some ad hoc experiments to produce configurations of complex lines with large 
packing diameters. For each pair {d,N), we used the Rankin bound to determine a lower limit on 
the feasibility parameter /x. Starting matrices were constructed with Algorithm S] using values of 
r ranging between 0.9 and 1.0. (Algorithm [J] typically fails for smaller values of r.) For values of 
the feasibility parameter between the minimal value and twice the minimal value, we performed 
5000 iterations of Algorithm [H and we recorded the largest packing diameter attained during these 
trials. 

Table [2] compares our results against the Rankin bound. We see that many of the complex line 
configurations have packing diameters much smaller than the Rankin bound, which is not surprising 
because the bound is usually not attainable. Some of our configurations fall within a thousandth 
of a degree of the bound, which is essentially optimal. 

Table [2] contains a few oddities. In P'^(C), the best packing diameter computed for = 
18, 19, . . . , 24 is worse than the packing diameter for A^ = 25. This configuration of 25 lines is 
an equiangular tight frame, which means that it is an optimal packing [TDJSOSt Table 1] . It seems 
likely that the optimal configurations for the preceding values of A^ are just subsets of the optimal 
arrangement of 25 lines. As before, it may be difficult to calculate this type of degenerate packing. 
A similar event occurs less dramatically at the pair (d, N) = (4, 13) and at the pairs (4, 17) and 
(4,18). 

Figure [T] compares the quality of the best real projective packings from |Slo04a| with the best 
complex projective packings that we obtained. It is natural that the complex packings are better 
than the real packings because the real projective space can be embedded isometrically into the 
complex projective space. But it is remarkable how badly the real packings compare with the 
complex packings. The only cases where the real and complex ensembles have the same packing 
diameter occur when the real configuration meets the Rankin bound. 

5.2. The Chordal Distance. Emboldened by this success with projective packings, we move on 
to packings of subspaces with respect to the chordal distance. Once again, we are able to use 
Sloane's tables for guidance in the real case. In the complex case, we fall back on the Rankin 
bound. 

For each triple {d,K,N), we determined a value for the feasibility parameter fi from the best 
packing diameter Sloane recorded for A^ subspaces in G{K, M'^), along with equation (j3.2p . We con- 
structed starting points using the modified version of Algorithm H] with r = ^/K, which represents 
no constraint. (We found that the alternating projection performed no better with initial config- 
urations generated from smaller values of r.) Then we executed Algorithm [1] with the calculated 
value of fj, for 5000 iterations. 
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Table [3] demonstrates how the best packings we obtained compare with Sloane's best packings. 
Many of our real configurations attained a squared packing diameter within 10"'^ of the best value 
Sloane recorded. Our algorithm was especially successful for smaller numbers of subspaces, but its 
performance began to flag as the number of subspaces approached 20. 

Table [3] contains several anomalies. For example, our configurations of = 11, 12, . . . , 16 sub- 
spaces in M'* yield worse packing diameters than the configuration of 17 subspaces. It turns out that 
this configuration of 17 subspaces is optimal, and Sloane's data show that the (putative) optimal 
arrangements of 11 to 16 subspaces are all subsets of this configuration. This is the same problem 
that occurred in some of our earlier experiments, and it suggests again that our algorithm has 
difficulty locating these degenerate configurations precisely. 

The literature contains very few experimental results on packing in complex Grassmannian man- 
ifolds equipped with chordal distance. To our knowledge, the only numerical work appears in two 
short tables from [ARUOlj . Therefore, we found it valuable to compare our results against the 
Rankin bound for subspace packings, which is derived in Section [H For reference, this bound 
requires that every configuration ^ of N subspaces in G{K, F'^) satisfy the inequality 

2 K{d-K) N 
pack(,ijord(^ ) ^ ■ 



d N -I 

This bound cannot always be met. In particular, the bound is attainable in the complex setting 
only if < d^. In the real setting, the bound requires that N < ^d{d + 1). When the bound is 
attained, each pair of subspaces in ^ is equidistant. 

We performed some ad hoc experiments to construct a table of packings in G{K, C^) equipped 
with the chordal distance. For each triple {d,K,N), we constructed random starting points using 
Algorithm m with r = ^/K (which represents no constraint). Then we used the Rankin bound 
to calculate a lower limit on the feasibility parameter fi. For this value of fi, we executed the 
alternating projection. Algorithm [H for 5000 iterations. 

The best packing diameters we obtained are listed in Table HI We see that there is a remarkable 
correspondence between the squared packing diameters of our configurations and the Rankin bound. 
Indeed, many of our packings are within 10""^ of the bound, which means that these configurations 
are essentially optimal. The algorithm was less successful as approached d^, which is an upper 
bound on the number N of subspaces for which the Rankin bound is attainable. 

Figure [2] compares the packing diameters of the best configurations in real and complex Grass- 
mannian spaces equipped with chordal distance. It is remarkable that both real and complex 
packings almost meet the Rankin bound for all N where it is attainable. Notice how the real pack- 
ing diameters fall off as soon as N exceeds ^d[d + 1). In theory, a complex configuration should 
always attain a better packing diameter than the corresponding real configuration because the real 
Grassmannian space can be embedded isometrically into the complex Grassmannian space. The 
figure shows that our best arrangements of 17 and 18 subspaces in G(2, C^) are actually slightly 
worse than the real arrangements calculated by Sloane. This indicates a failure of the alternating 
projection algorithm. 

5.3. The Spectral Distance. Next, we consider how to compute Grassmannian packings with 
respect to the spectral distance. This investigation requires some small modifications to the al- 
gorithm, which are described in the next subsection. Afterward, we provide the results of some 
numerical experiments. 

5.3.1. Modifications to Algorithm. To construct packings with respect to the spectral distance, we 
tread a familiar path. Suppose that we wish to produce a configuration of subspaces in G{K, C^) 
with a packing diameter p. The feasibility problem requires that 

max ||X*,X„||2 2 < (5-1) 
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where fi = \J\ — p^. This leads to the convex structural constraint set 

= {H G C^'^^^^ :H = H\ Hnn = I forn = 1, 2, . . . , iV, and 

||-?^mn|l2 2 — ^ 

The spectral constraint set is the same as before. The next proposition shows how to find the matrix 
in ^ closest to an initial matrix. In preparation, define the truncation operator [x]^ = min{x, /x} 
for numbers, and extend it to matrices by applying it to each component. 

Proposition 9. Let G be an Hermitian matrix. With respect to the Frobenius norm, the unique 
matrix in J^{fi) nearest to G has a block identity diagonal. If the off-diagonal block Gmn has a 
singular value decomposition UmnCmnVmnj then 

GjYin if II G^mn II 2 2 — 

Umn [Cmn]^l otherwise. 

Proof. To determine the {m, n) off-diagonal block of the solution matrix H, we must solve the 
optimization problem 

111 1 1 2 

min^ 2 11^ — Gmnllp subject to ||A||2 2<Ai- 

The Frobenius norm is strictly convex and the spectral norm is convex, so this problem has a unique 
solution. 

Let <t(-) return the vector of decreasingly ordered singular values of a matrix. Suppose that G^n 
has the singular value decomposition Gmn = U{diag a{Gmn)}V* ■ The constraint in the optimiza- 
tion problem depends only on the singular values of A, and so the Hoffman-Wielandt Theorem for 
singular values |HJ85j allows us to check that the solution has the form A = U {diag cr {A)}V * . 

To determine the singular values ^ = cr{A) of the solution, we must solve the (strictly) convex 
program 

min^ i 11^ - cr(Gmn)|l2 subject to < fJ- 
An easy application of Karush-Kuhn-Tucker theory [Roc70] proves that the solution is obtained 
by truncating the singular values of Gmn that exceed fj,. □ 

5.3.2. Numerical Results. To our knowledge, there are no numerical studies of packing in Grass- 
mannian spaces equipped with spectral distance. To gauge the quality of our results, we compare 
them against the upper bound of Corollary [71 In the real or complex setting, a configuration J^T of 
subspaces in G{K, ¥'^) with respect to the spectral distance must satisfy the bound 

, d- K N 

Pack,pcc(^) < ^ WZTJ- 

In the real case, the bound is attainable only if A < ^d{d+l) — ^ K [K -|- 1) -|- 1, while attainment 
in the complex case requires that N < d'^ — + 1 |LS73j . When a configuration meets the bound, 
the subspaces are not only equidistant but also equi-isoclinic. That is, all principal angles between 
all pairs of subspaces are identical. 

We performed some limited ad hoc experiments in an effort to produce good configurations of 
subspaces with respect to the spectral distance. We constructed random starting points using the 
modified version of Algorithm [J] with r = 1, which represents no constraint. (Again, we did not 
find that smaller values of r improved the performance of the alternating projection.) Prom the 
Rankin bound, we calculated the smallest possible value of the feasibility parameter fi. For values 
of /X ranging from the minimal value to twice the minimal value, we ran the alternating projection. 
Algorithm [H for 5000 iterations, and we recorded the best packing diameters that we obtained. 

Table [5] displays the results of our calculations. We see that some of our configurations essentially 
meet the Rankin Bound, which means that they are equi-isoclinic. It is clear that alternating 
projection also succeeds reasonably well for this packing problem. 
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The most notable pathology in the table occurs for configurations of 8 and 9 subspaces in 
G(3, M^). In these cases, the algorithm always yielded arrangements of subspaces with a zero 
packing diameter, which implies that two of the subspaces intersect nontrivially. Nevertheless, we 
were able to construct random starting points with a nonzero packing diameter, which means that 
the algorithm is making the initial configuration worse. We do not understand the reason for this 
failure. 

Figure [3] makes a graphical comparison between the real and complex subspace packings. On 
the whole, the complex packings are much better than the real packings. For example, every 
configuration of subspaces in G(2,C^) nearly meets the Rankin bound, while just two of the real 
configurations achieve the same distinction. In comparison, it is curious how few arrangements in 
G(2,C^) come anywhere near the Rankin bound. 

5.4. The Fubini— Study Distance. When we approach the problem of packing in Grassmannian 
manifolds equipped with the Fubini-Study distance, we are truly out in the wilderness. To our 
knowledge, the literature contains neither experimental nor theoretical treatments of this ques- 
tion. Moreover, we are not presently aware of general upper bounds on the Fubini-Study packing 
diameter that we might use to assay the quality of a configuration of subspaces. Nevertheless, 
we attempted a few basic experiments. The investigation entails some more modifications to the 
algorithm, which are described below. Afterward, we go over our experimental results. We view 
this work as very preliminary. 

5.4.1. Modifications to Algorithm. Suppose that we wish to construct a configuration of sub- 
spaces whose Fubini-Study packing diameter exceeds p. The feasibility condition is 

max |detXj^X.„| < /_i (5.2) 
where /u = cos p. This leads to the structural constraint set 

^(/i) = {H G C^^x^^ :H = H\ = I for n = 1, 2, . . . , iV, and 

I det Hmn I < for all m ^ n}. 

Unhappily, this set is no longer convex. To produce a nearest matrix in J^, we must solve a 
nonlinear programming problem. The following proposition describes a numerically favorable for- 
mulation. 

Proposition 10. Let G he an Hermitian matrix. Suppose that the off-diagonal block Gmn has 
singular value decomposition UmnCmnVmn- -^^i Cmn = diagCmn, and find a (real) vector Xmn that 
solves the optimization problem 

min i ||exp(a;) — Cmnllo subject to e* x < log p. 

X 

In Frobenius norm, a matrix H from J^{p) that is closest to G has a block-identity diagonal and 
off-diagonal blocks 

H =[ ^"'"^ i/|detGmn| < P, and 

™" 1 J7mn{diag(expa;„„)}V„*„ otherwise. 

We use exp(-) to denote the componentwise exponential of a vector. One may establish that the 
optimization problem is not convex by calculating the Hessian of the objective function. 

Proof. To determine the (m, n) off-diagonal block of the solution matrix we must solve the 
optimization problem 

111 1 1 2 

miiiyi 2 11^ — Gmnllp subject to |detA|</x. 
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We may reformulate this problem as 

miiiA ^ ||A - GmnWl subject to ^^^^ log o-fc(A) < log^. 

A familiar argument proves that the solution matrix has the same left and right singular vectors 
as Gmn- To obtain the singular values ^ = o'(A) of the solution, we consider the mathematical 
program 

min^ i 11^ - <j{Gm.n)\\2 subject to X^^^-^log^fe < log//. 
Change variables to complete the proof. □ 

5.4.2. Numerical Experiments. We implemented the modified version of Algorithm [1] in Matlab, 
using the built-in nonlinear programming software to solve the optimization problem required by 
the proposition. For a few triples (d, K, N), we ran 100 to 500 iterations of the algorithm for various 
values of the feasibility parameter ^. (Given the exploratory nature of these experiments, we found 
that the implementation was too slow to increase the number of iterations.) 

The results appear in Table [6j For small values of N, we find that the packings exhibit the 
maximum possible packing diameter tt/2, which shows that the algorithm is succeeding in these 
cases. For larger values of N, we are unable to judge how close the packings might decline from 
optimal. 

Figure H] compares the quality of our real packings against our complex packings. In each case, 
the complex packing is at least as good as the real packing, as we would expect. The smooth 
decline in the quality of the complex packings suggests that there is some underlying order to the 
packing diameters, but it remains to be discovered. 

To perform large-scale experiments, it will probably be necessary to tailor an algorithm that can 
solve the nonlinear programming problems more quickly. It may also be essential to implement 
the alternating projection in a programming environment more efficient than Matlab. Therefore, a 
detailed study of packing with respect to the Fubini-Study distance must remain a topic for future 
research. 

6. Discussion 

6.1. Subspace Packing in Wireless Communications. Configurations of subspaces arise in 
several aspects of wireless communication, especially in systems with multiple transmit and receive 
antennas. The intuition behind this connection is that the transmitted and received signals in a 
multiple antenna system are connected by a matrix transformation, or matrix channel. 

Subspace packings occur in two wireless applications: noncoherent communication and in sub- 
space quantization. The noncoherent application is primarily of theoretical interest, while subspace 
quantization has a strong impact on practical wireless systems. Grassmannian packings appear in 
these situations due to an assumption that the matrix channel should be modeled as a complex 
Gaussian random matrix. 

In the noncoherent communication problem, it has been shown that, from an information- 
theoretic perspective, under certain assumptions about the channel matrix, the optimum transmit 
signal corresponds to a packing in G{K, C^) where K corresponds to the minimum of the number of 
transmit and receive antennas and d corresponds to the number of consecutive samples over which 
the channel is constant |ZT021 IHMOOj . In other words, the number of subspaces K is determined by 
the system configuration, while d is determined by the carrier frequency and the degree of mobility 
in the propagation channel. 

On account of this application, several papers have investigated the problem of finding packings 
in Grassmannian manifolds. One approach for the case oi K = 1 is presented in [HMOOj . This paper 
proposes a numerical algorithm for finding line packings, but it does not discuss its properties or 
connect it with the general subspace packing problem. Another approach, based on discrete Fourier 
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transform matrices, appears in [HMR^OO]. This construction is both structured and flexible, but 
it does not lead to optimal packings. The paper |ARU01j studies Grassmannian packings in detail, 
and it contains an algorithm for finding packings in the complex Grassmannian manifold equipped 
with chordal distance. This algorithm is quite complex: it uses surrogate functionals to solve a 
sequence of relaxed nonlinear programs. The authors tabulate several excellent chordal packings, 
but it is not clear whether their method generalizes to other metrics. 

The subspace quantization problem also leads to Grassmannian packings. In multiple-antenna 
wireless systems, one must quantize the dominant subspace in the matrix communication channel. 
Optimal quantizers can be viewed as packings in G{K, C^), where K is the dimension of the subspace 
and d is the number of transmit antennas. The chordal distance, the spectral distance, and the 
Fubini-Study distance are all useful in this connection [LHJ051 ILJOSj . This literature does not 
describe any new algorithms for constructing packings; it leverages results from the noncoherent 
communication literature. Communication strategies based on quantization via subspace packings 
have been incorporated into at least one recent standard |Wir05] . 

6.2. Conclusions. We have shown that the alternating projection algorithm can be used to solve 
many different packing problems. The method is easy to understand and to implement, even while 
it is versatile and powerful. In cases where experiments have been performed, we have often been 
able to match the best packings known. Moreover, we have extended the method to solve problems 
that have not been studied numerically. Using the Rankin bounds, we have been able to show 
that many of our packings are essentially optimal. It seems clear that alternating projection is an 
effective numerical algorithm for packing. 

Appendix A. Tammes' Problem 

The alternating projection method can also be used to study Tammes' Problem of packing points 
on a sphere |Tam30j . This question has received an enormous amount of attention over the last 
75 years, and extensive tables of putatively optimal packings are available |Slo04b] . This appendix 
offers a brief treatment of our work on this problem. 

A.l. Modifications to Algorithm. Suppose that we wish to produce a configuration of N points 
on the unit sphere §>'^~^ with a packing diameter p. The feasibility problem requires that 

max: {x^, Xn) < fJ. (A.l) 

where p = \/l — p^. This leads to the convex structural constraint set 

J^ifi) = {H £ R'^'"^ :H = H*, hnn = lior n = 1,2,..., N, and 

— 1 < ^ fJ- for all m ^ n}. 

The spectral constraint set is the same as before. The associated matrix nearness problem is trivial 
to solve. 

Proposition 11. Let G be a real, symmetric matrix. With respect to Frobenius norm, the unique 
matrix in J^{n) closest to G has a unit diagonal and off-diagonal entries that satisfy 

1 ) 9mn ^ 1 5 

9mnt 1 !^ 9Tnn ^ l^i and 

^, II- < 9mn- 
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A. 2. Numerical Results. Tammes' Problem has been studied for 75 years, and many putatively 
optimal configurations are available. Therefore, we attempted to produce packings whose maximum 
inner product /i fell within 10"^ of the best value tabulated by N. J. A. Sloane and his colleagues 
|Slo04bj . This resource draws from all the experimental and theoretical work on Tammes' Problem, 
and it should be considered the gold standard. 

Our experimental setup echoes the setup for real projective packings. We implemented the 
algorithms in Matlab, and we performed the following experiment for pairs (d, A^) with d = 3, 4, 5 
and A = 4, 5, ... , 25. First, we computed the putatively optimal maximum inner product fi using 
the data from |Slo04b| . In each of 10 trials, we constructed a starting matrix using Algorithm [J] with 
parameters r = 0.9 and T = 10, 000. Then, we executed the alternating projection, Algorithm [H 
with the calculated value of and the maximum number of iterations set to T = 5000. We stopped 
the alternating projection in Step 4 if the iterate G^*-* contained no off-diagonal entry greater than 
fi + 10~^ and proceeded with Step 6. After 10 trials, we recorded the largest packing diameter 
attained, as well as the average value of the packing diameter. We also recorded the average 
number of iterations the alternating projection required during each trial. 

Table [7] provides the results of this experiment. The most striking feature of Table [7] is that 
the best configurations returned by alternating projection consistently attain packing diameters 
that fall hundredths or thousandths of a degree away from the best packing diameters recorded by 
Sloane. If we examine the maximum inner product in the configuration instead, the difference is 
usually on the order of 10~^ or 10~^, which we expect based on our stopping criterion. The average- 
case results are somewhat worse. Nevertheless, the average configuration returned by alternating 
projection typically attains a packing diameter only several tenths of a degree away from optimal. 

A second observation is that the alternating projection algorithm typically performs better when 
the number of points N is small. The largest errors are all clustered at larger values of N. A 
corollary observation is that the average number of iterations per trial tends to increase with the 
number of points. We believe that the explanation for these phenomena is that Tammes' Problem 
has a combinatorial regime, where solutions have a lot of symmetry and structure, and a random 
regime, where the solutions have very little order. The algorithm typically seems to perform better 
in the combinatorial regime, although it fails for certain unusually structured ensembles. 

This claim is supported somewhat by theoretical results for d = 3. Optimal configurations have 
only been established for N = 1, 2, . . . , 12 and N = 24. Of these, the cases N = 1, 2, 3 are trivial. 
The cases A^ = 4, 6, 8, 12, 24 fall from the vertices of various well-known polyhedra. The cases 
A^ = 5,11 are degenerate, obtained by leaving a point out of the solutions for A = 6,12. The 
remaining cases involve complicated constructions based on graphs |EZ01j . The algorithm was 
able to calculate the known optimal configurations to a high order of accuracy, but it generally 
performed slightly better for the nondegenerate cases. 

On the other hand, there is at least one case where the algorithm failed to match the optimal 
packing diameter, even though the optimal configuration is highly symmetric. The best arrange- 
ment of 24 points on locates them at vertices of a very special polytope called the 24-cell |Slo04bj . 
The best configuration produced by the algorithm has a packing diameter 1.79° worse. It seems 
that this optimal configuration is very difficult for the algorithm to find. Less dramatic failures 
occurred at pairs {d,N) = (3,25), (4,14), (4,25), (5,22), and (5,23). But in each of these cases, 
our best packing declined more than a tenth of a degree from the best recorded. 
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Appendix B. Tables and Figures 

Our experiments resulted in tables of packing diameters. We did not store the configurations 
produced by the algorithm. The Matlab code that produced these data is available on request from 
jtroppOacm. caltech.edu. 

These tables and figures are intended only to describe the results of our experiments; it is likely 
that many of the packing diameters could be improved with additional effort. In all cases, we 
present the results of calculations for the stated problem, even if we obtained a better packing 
by solving a different problem. For example, a complex packing should always improve on the 
corresponding real packing. If the numbers indicate otherwise, it just means that the complex 
experiment yielded an inferior result. As a second example, the optimal packing diameter must not 
decrease as the number of points increases. When the numbers indicate otherwise, it means that 
running the algorithm with more points yielded a better result than running it with fewer. These 
failures may reflect the difficulty of various packing problems. 
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Table 1. Packing in real projective spaces: For collections of iV points in the 
real projective space P'^~^(]R), this table lists the best packing diameter (in degrees) 
and the average packing diameter (in degrees) obtained during ten random trials of 
the alternating projection algorithm. The error columns record how far our results 
decline from the putative optimal packings (NJAS) reported in |Slo04aj . The last 
column gives the average number of iterations of alternating projection per trial 
before the termination condition is met. 
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. . . continued 





Packing diameters (Degrees) 
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Table 2. Packing in complex projective spaces: We compare our best con- 
figurations (DHST) of A'^ points in the complex projective space F'^~^(C) against 
the Rankin bound (j4.ip . The packing diameter of an ensemble is measured as the 
acute angle (in degrees) between the closest pair of lines. The final column shows 
how far our configurations fall short of the bound. 
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. . . continued 





Packing diameters (Degrees) 
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Figure 1. Real and Complex Projective Packings: These three graphs com- 
pare the packing diameters attained by configm'ations in real and complex projective 
spaces with d = 3,4, 5. The circles indicate the best real packings obtained by Sloane 
and his colleagues |Slo04a| . The crosses indicate the best complex packings produced 
by the authors. Rankin's upper bound (j4.ip is depicted in gray. The dashed vertical 
line marks the largest number of real lines for which the Rankin bound is attainable, 
while the solid vertical line marks the maximum number of complex lines for which 
the Rankin bound is attainable. 
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Table 3. Packing in real Grassmannians with chordal distance: We com- 
pare our best configurations (DHST) of N points in G{K, M'^) against the best 
packings (NJAS) reported in |Slo04a] . The squared packing diameter is the squared 
chordal distance ()2.ip between the closest pair of subspaces. The last column lists 
the difference between the columns (NJAS) and (DHST). 
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Table 4. Packing in complex Grassmannians with chordal distance: We 
compare our best configurations (DHST) of N points in G{K, C"') against tlie Rankin 
bound, equation (14. ip . The squared packing diameter is calculated as the squared 
chordal distance ()2.ip between the closest pair of subspaces. The final column shows 
how much the computed ensemble declines from the Rankin bound. When the bound 
is met, all pairs of subspaces are equidistant. 
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. . . continued 
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. . . continued 
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Figure 2. Packing in Grassmannians with chordal distance: This figure 
shows the packing diameters of N points in the Grassmannian G{K, W^) equipped 
with the chordal distance. The circles indicate the best real packings (F = M) 
obtained by Sloane and his colleagues |Slo04aj . The crosses indicate the best complex 
packings (F = C) produced by the authors. Rankin's upper bound (|4.ip appears in 
gray. The dashed vertical line marks the largest number of real subspaces for which 
the Rankin bound is attainable, while the solid vertical line marks the maximum 
number of complex subspaces for which the Rankin bound is attainable. 
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Table 5. Packing in Grassmannians with spectral distance: We compare 
our best real (F = M) and complex (F = C) packings in G{K,¥'^) against the 
Rankin bound, equation ()4.2p . The squared packing diameter of a configuration is 
the squared spectral distance ()2.2p between the closest pair of subspaces. When the 
Rankin bound is met, all pairs of subspaces are equi-isoclinic. The algorithm failed 
to produce any configurations of 8 or 9 subspaces in G(3, M^) with nontrivial packing 
diameters. 
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Figure 3. Packing in Grassmannians with spectral distance: This figure 
shows the packing diameters of N points in the Grassmannian G{K, W^) equipped 
with the spectral distance. The circles indicate the best real packings (F = M) 
obtained by the authors, while the crosses indicate the best complex packings (F = 
C) obtained. The Rankin bound (j4.2p is depicted in gray. The dashed vertical line 
marks an upper bound on largest number of real subspaces for which the Rankin 
bound is attainable according to Theorem [51 
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Table 6. Packing in Grassmannians with Fubini-Study distance: Our best 
real packings (F = R) compared with our best complex packings (F = C) in the space 
G{K, F*^) . The packing diameter of a configuration is the Fubini-Study distance (|2.3p 
between the closest pair of subspaces. Note that we have scaled the distance by 2/7r 
so that it ranges between zero and one. 
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Figure 4. Packing in Grassmannians with Fubini-Study distance: This 
figure shows the packing diameters of points in the Grassmannian G{K,¥'^) 
equipped with the Fubini-Study distance. The circles indicate the best real pack- 
ings (F = M) obtained by the authors, while the crosses indicate the best complex 
packings (F = C) obtained. 
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Table 7. Packing on spheres: For collections of N points on the (d — 1)- 
dimensional sphere, this table lists the best packing diameter and the average 
packing diameter obtained during ten random trials of the alternating projection 
algorithm. The error columns record how far our results decline from the putative 
optimal packings (NJAS) reported in [Slo04b| . The last column gives the average 
number of iterations of alternating projection per trial. 
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. . . continued 
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